import sympy as sym
import sympy as sym
import numpy as np
from sympy.physics.mechanics import dynamicsymbols, init_vprinting
from sympy import integrate, symbols, Matrix, Eq, solve, simplify, cos, sin, det, pi, lambdify, BlockMatrix
init_vprinting()
def skew(w_):
return Matrix([[0, -w_[2], w_[1]],
[w_[2], 0, -w_[0]],
[-w_[1], w_[0], 0]])
def unskew(w_):
return Matrix([w_[2,1], w_[0,2], w_[1,0]])
def hat(V_):
v = V_[0:3,0]
w = V_[3:6,0]
bottomones = Matrix([[0,0,0,1]])
return skew(w).row_join(v).col_join(bottomones)
def unhat(Vhat_):
what = Vhat_[0:3,0:3]
v = Vhat_[0:3,3]
return v.col_join(unskew(what))
def Rot(th_):
return Matrix([[cos(th_), -sin(th_), 0],
[sin(th_), cos(th_), 0],
[0,0,1]])
def hom(v_):
return v_.col_join(Matrix([1]))
def unhom(v_):
return Matrix([v_[0,:],v_[1,:],v_[2,:]])
def g(x_,y_,th_):
translation = Matrix([x_,y_,0]) # vector p
bottomones = Matrix([[0,0,0,1]]) # bottom row of g
return Rot(th_).row_join(translation).col_join(bottomones)
def ginv(g_):
R = g_[0:3,0:3]
p = g_[0:3,3]
bottomones = Matrix([[0,0,0,1]])
return (R.T).row_join(-R.T*p).col_join(bottomones)
def Itens(m_,I_):
mat = m_*sym.eye(6)
mat[3:6,3:6] = I_
return mat
xb,yb,thb,xj,yj,thj = dynamicsymbols(r'x_{box} y_{box} \theta_{box} x_{jack} y_{jack} \theta_{jack}')
M,m,L,l,grav,t = symbols(r'M m L l g t')
q = Matrix([xb,yb,thb,xj,yj,thj])
qdot = q.diff(t)
qddot = qdot.diff(t)
lam = symbols(r'\lambda')
givens = {M:10,m:1,L:1,l:.2,grav:9.8}
gwa = g(xb,yb,0)
gab = g(0,0,thb)
gwb = gwa*gab # TRANSFORMATION FOR CENTER OF BOX
Vwb = simplify(unhat(ginv(gwb)*gwb.diff(t)))
rb = sym.Matrix([0,0,0])
gwc = g(xj,yj,0)
gcd = g(0,0,thj)
gwd = gwc*gcd # TRANSFORMATION FOR CENTER OF JACK
Vwd = simplify(unhat(ginv(gwd)*gwd.diff(t)))
rd = sym.Matrix([0,0,0])
gbd = simplify(ginv(gwb)*gwd) # TRANSFORM JACK FRAME TO BOX FRAME
J = (1/6)*M*(L**2)
Ib = sym.Matrix([[0, 0, 0],
[0, 0, 0],
[0, 0, J]])
j = (1/6)*m*1**2
Id = sym.Matrix([[0, 0, 0],
[0, 0, 0],
[0, 0, j]])
KE = simplify(((1/2)*Vwb.T*Itens(M,Ib)*Vwb + (1/2)*Vwd.T*Itens(m,Id)*Vwd)[0])
V = simplify((M*grav*Matrix([0,1,0,0]).T*gwb*hom(rb) + m*grav*Matrix([0,1,0,0]).T*gwc*hom(rd))[0])
thbd = pi/20 + (pi/45)*sin(t/20)**2
Fthb = -20*(thb-thbd)
F = Matrix([0,1.1 * M*grav,Fthb,0,0,0]).subs(givens)
Lag = (KE-V).subs(givens)
ELlhs = Lag.diff(qdot).diff(t) - Lag.diff(q)
ELeq = Eq(ELlhs, F)
soln_reg = solve(ELeq,qddot)
soln_reg
xbddot_fn = lambdify([[q[0],q[1],q[2],q[3],q[4],q[5],qdot[0],qdot[1],qdot[2],qdot[3],qdot[4],qdot[5]],t],soln_reg[qddot[0]])
ybddot_fn = lambdify([[q[0],q[1],q[2],q[3],q[4],q[5],qdot[0],qdot[1],qdot[2],qdot[3],qdot[4],qdot[5]],t],soln_reg[qddot[1]])
thbddot_fn = lambdify([[q[0],q[1],q[2],q[3],q[4],q[5],qdot[0],qdot[1],qdot[2],qdot[3],qdot[4],qdot[5]],t],soln_reg[qddot[2]])
xjddot_fn = lambdify([[q[0],q[1],q[2],q[3],q[4],q[5],qdot[0],qdot[1],qdot[2],qdot[3],qdot[4],qdot[5]],t],soln_reg[qddot[3]])
yjddot_fn = lambdify([[q[0],q[1],q[2],q[3],q[4],q[5],qdot[0],qdot[1],qdot[2],qdot[3],qdot[4],qdot[5]],t],soln_reg[qddot[4]])
thjddot_fn = lambdify([[q[0],q[1],q[2],q[3],q[4],q[5],qdot[0],qdot[1],qdot[2],qdot[3],qdot[4],qdot[5]],t],soln_reg[qddot[5]])
lag_fn = [xbddot_fn,ybddot_fn,thbddot_fn,xjddot_fn,yjddot_fn,thjddot_fn]
# Box
A = L/2
B = L/2
C = -L/2
D = -L/2
# Jack
rd0 = sym.Matrix([l/2,0,0])
rd1 = sym.Matrix([0,l/2,0])
rd2 = sym.Matrix([-l/2,0,0])
rd3 = sym.Matrix([0,-l/2,0])
phi0A = A - (Matrix([1,0,0,0]).T * gbd * hom(rd0))[0]
phi0B = B - (Matrix([0,1,0,0]).T * gbd * hom(rd0))[0]
phi0C = (Matrix([1,0,0,0]).T * gbd * hom(rd0))[0] - C
phi0D = (Matrix([0,1,0,0]).T * gbd * hom(rd0))[0] - D
phi1A = A - (Matrix([1,0,0,0]).T * gbd * hom(rd1))[0]
phi1B = B - (Matrix([0,1,0,0]).T * gbd * hom(rd1))[0]
phi1C = (Matrix([1,0,0,0]).T * gbd * hom(rd1))[0] - C
phi1D = (Matrix([0,1,0,0]).T * gbd * hom(rd1))[0] - D
phi2A = A - (Matrix([1,0,0,0]).T * gbd * hom(rd2))[0]
phi2B = B - (Matrix([0,1,0,0]).T * gbd * hom(rd2))[0]
phi2C = (Matrix([1,0,0,0]).T * gbd * hom(rd2))[0] - C
phi2D = (Matrix([0,1,0,0]).T * gbd * hom(rd2))[0] - D
phi3A = A - (Matrix([1,0,0,0]).T * gbd * hom(rd3))[0]
phi3B = B - (Matrix([0,1,0,0]).T * gbd * hom(rd3))[0]
phi3C = (Matrix([1,0,0,0]).T * gbd * hom(rd3))[0] - C
phi3D = (Matrix([0,1,0,0]).T * gbd * hom(rd3))[0] - D
delphi0A = phi0A.diff(q)
delphi0B = phi0B.diff(q)
delphi0C = phi0C.diff(q)
delphi0D = phi0D.diff(q)
delphi1A = phi1A.diff(q)
delphi1B = phi1B.diff(q)
delphi1C = phi1C.diff(q)
delphi1D = phi1D.diff(q)
delphi2A = phi2A.diff(q)
delphi2B = phi2B.diff(q)
delphi2C = phi2C.diff(q)
delphi2D = phi2D.diff(q)
delphi3A = phi3A.diff(q)
delphi3B = phi3B.diff(q)
delphi3C = phi3C.diff(q)
delphi3D = phi3D.diff(q)
xb_t = sym.symbols(r'x_{box}(\tau^-)')
yb_t = sym.symbols(r'y_{box}(\tau^-)')
thb_t = sym.symbols(r'\theta_{box}(\tau^-)')
xj_t = sym.symbols(r'x_{jack}(\tau^-)')
yj_t = sym.symbols(r'y_{jack}(\tau^-)')
thj_t = sym.symbols(r'\theta_{jack}(\tau^-)')
xbdot_t = sym.symbols(r'\dot{x}_{box}(\tau^-)')
ybdot_t = sym.symbols(r'\dot{y}_{box}(\tau^-)')
thbdot_t = sym.symbols(r'\dot{\theta}_{box}(\tau^-)')
xjdot_t = sym.symbols(r'\dot{x}_{jack}(\tau^-)')
yjdot_t = sym.symbols(r'\dot{y}_{jack}(\tau^-)')
thjdot_t = sym.symbols(r'\dot{\theta}_{jack}(\tau^-)')
q_t = Matrix([xb_t,yb_t,thb_t,xj_t,yj_t,thj_t])
qdot_t = Matrix([xbdot_t, ybdot_t, thbdot_t, xjdot_t, yjdot_t, thjdot_t])
qsubs = {}
for i in range(len(q)):
qsubs[q[i]] = q_t[i]
qsubs[qdot[i]] = qdot_t[i]
# need p_t, delphi_t, H_t
dLdqdot = Lag.diff(qdot)
p_t = sym.nsimplify(dLdqdot.subs(qsubs))
print('p_t =')
display(p_t)
delphi0A_t = delphi0A.subs(qsubs)
delphi0B_t = delphi0B.subs(qsubs)
delphi0C_t = delphi0C.subs(qsubs)
delphi0D_t = delphi0D.subs(qsubs)
delphi1A_t = delphi1A.subs(qsubs)
delphi1B_t = delphi1B.subs(qsubs)
delphi1C_t = delphi1C.subs(qsubs)
delphi1D_t = delphi1D.subs(qsubs)
delphi2A_t = delphi2A.subs(qsubs)
delphi2B_t = delphi2B.subs(qsubs)
delphi2C_t = delphi2C.subs(qsubs)
delphi2D_t = delphi2D.subs(qsubs)
delphi3A_t = delphi3A.subs(qsubs)
delphi3B_t = delphi3B.subs(qsubs)
delphi3C_t = delphi3C.subs(qsubs)
delphi3D_t = delphi3D.subs(qsubs)
# print('delphi_t=')
# delphi_t = delphi.subs(qsubs)
# display(delphi_t)
H_t = ((dLdqdot.T*qdot)[0] - Lag).subs(qsubs)
print('Hamiltonian =')
display(H_t)
p_t =
Hamiltonian =
xbdot_p = sym.symbols(r'\dot{x}_{box}(\tau^+)')
ybdot_p = sym.symbols(r'\dot{y}_{box}(\tau^+)')
thbdot_p = sym.symbols(r'\dot{\theta}_{box}(\tau^+)')
xjdot_p = sym.symbols(r'\dot{x}_{jack}(\tau^+)')
yjdot_p = sym.symbols(r'\dot{y}_{jack}(\tau^+)')
thjdot_p = sym.symbols(r'\dot{\theta}_{jack}(\tau^+)')
qdot_p = Matrix([xbdot_p,ybdot_p,thbdot_p,xjdot_p,yjdot_p,thjdot_p])
qdotpsubs = {qdot_t[i]:qdot_p[i] for i in range(len(qdot_t))}
eq0A = simplify(Eq(p_t.subs(qdotpsubs) - p_t, lam*delphi0A_t))
eq0B = simplify(Eq(p_t.subs(qdotpsubs) - p_t, lam*delphi0B_t))
eq0C = simplify(Eq(p_t.subs(qdotpsubs) - p_t, lam*delphi0C_t))
eq0D = simplify(Eq(p_t.subs(qdotpsubs) - p_t, lam*delphi0D_t))
eq1A = simplify(Eq(p_t.subs(qdotpsubs) - p_t, lam*delphi1A_t))
eq1B = simplify(Eq(p_t.subs(qdotpsubs) - p_t, lam*delphi1B_t))
eq1C = simplify(Eq(p_t.subs(qdotpsubs) - p_t, lam*delphi1C_t))
eq1D = simplify(Eq(p_t.subs(qdotpsubs) - p_t, lam*delphi1D_t))
eq2A = simplify(Eq(p_t.subs(qdotpsubs) - p_t, lam*delphi2A_t))
eq2B = simplify(Eq(p_t.subs(qdotpsubs) - p_t, lam*delphi2B_t))
eq2C = simplify(Eq(p_t.subs(qdotpsubs) - p_t, lam*delphi2C_t))
eq2D = simplify(Eq(p_t.subs(qdotpsubs) - p_t, lam*delphi2D_t))
eq3A = simplify(Eq(p_t.subs(qdotpsubs) - p_t, lam*delphi3A_t))
eq3B = simplify(Eq(p_t.subs(qdotpsubs) - p_t, lam*delphi3B_t))
eq3C = simplify(Eq(p_t.subs(qdotpsubs) - p_t, lam*delphi3C_t))
eq3D = simplify(Eq(p_t.subs(qdotpsubs) - p_t, lam*delphi3D_t))
p_eq = [eq0A,eq0B,eq0C,eq0B,
eq1A,eq1B,eq1C,eq1B,
eq2A,eq2B,eq2C,eq2B,
eq3A,eq3B,eq3C,eq3B]
Heq = sym.simplify(Eq(H_t.subs(qdotpsubs) - H_t, 0))
print('Hamiltonian =')
display(Heq)
Hamiltonian =
def impact_update(s,p_eq):
subvar = {xb_t: s[0], yb_t: s[1], thb_t: s[2], xj_t:s[3], yj_t:s[4], thj_t:s[5],
xbdot_t: s[6], ybdot_t: s[7], thbdot_t: s[8], xjdot_t:s[9], yjdot_t:s[10], thjdot_t:s[11]}
p_eq_new = p_eq.subs(subvar).subs(givens)
Heq_new = Heq.subs(subvar)
imp_solns = solve((p_eq_new,Heq_new),(xbdot_p,ybdot_p,thbdot_p,xjdot_p,yjdot_p,thjdot_p,lam),dict=True)
finalsoln = [x for x in imp_solns if np.abs(x[lam]) >= 1e-5][0]
# display(imp_solns)
# print(finalsoln[lam])
sdot = np.array([finalsoln[qdot_p[i]] for i in range(len(qdot_p))]).T
return np.hstack([s[0:int(len(s)/2)],sdot])
print(impact_update(np.array([0,0,0,np.pi,0,0,-10,-1,-1,-1,-12,-1]),eq0A))
[0.0 0.0 0.0 3.141592653589793 0.0 0.0 -8.36363636363636 -1.00000000000000 -1.00000000000000 -17.3636363636364 -12.0000000000000 -1.00000000000000]
# DEBUGGING PLACES IN CODE WHERE SIMULATE GETS STUCK CONTINUOUSLY DETECTING IMPACT
test = np.array([ 0.00546604, 0.04357109, -0.0417529, -0.05466036, -0.43571085, 0.57955032,
-0.13760388, 0.23395788, -1.3991147, 1.37603882, -2.33957877, 0.11203142])
test2 = np.array([ 4.88632187e-04, 7.20929321e-02, 2.22022370e-01, -4.88632187e-03,
-4.20929321e-01, 1.90617991e+00, 1.08700449e-02, 3.23079384e-01,
8.93480300e-01, -1.08700449e-01, -3.23079384e+00, -2.17782093e+02])
impact_update(test2,p_eq[11])
array([0.000488632187, 0.0720929321, 0.22202237, -0.00488632187,
-0.420929321, 1.90617991, 0.242088135334012, -0.701170244417902,
1.53921976364961, -2.42088135334012, 7.01170244417902,
-217.069430120562], dtype=object)
def impact_condition(s, phi_fn):
return phi_fn(s[0:int(len(s)/2)]) < 0
# check all impacts are positive inside
center = np.array([0,0,0,0,0,0])
phis = [phi0A,phi0B,phi0C,phi0D,
phi1A,phi1B,phi1C,phi1D,
phi2A,phi2B,phi2C,phi2D,
phi3A,phi3B,phi3C,phi3D]
phi_fns = []
for p in phis:
phi_fns.append(lambdify([[q[0],q[1],q[2],q[3],q[4],q[5]]],p.subs(givens)))
for p in phi_fns:
# check that all phis are positive when the jack is inside the box
print(p(center))
0.4 0.5 0.6 0.5 0.5 0.4 0.5 0.6 0.6 0.5 0.4 0.5 0.5 0.6 0.5 0.4
# Hamiltonian
H = Lag.diff(qdot).dot(qdot) - Lag
ham_fn = lambdify([[q[0],q[1],q[2],q[3],q[4],q[5],qdot[0],qdot[1],qdot[2],qdot[3],qdot[4],qdot[5]]], H)
def integrate(f, xt, dt, t):
k1 = dt * f(xt,t)
k2 = dt * f(xt+k1/2.,t)
k3 = dt * f(xt+k2/2.,t)
k4 = dt * f(xt+k3,t)
new_xt = xt + (1/6.) * (k1+2.0*k2+2.0*k3+k4)
return new_xt
def sim_impact(f, x0, tspan, dt, integrate_f, impact_update_f, impact_cond, list_of_phi_fn, list_of_p_eq):
N = int((max(tspan)-min(tspan))/dt)
x = np.copy(x0)
tvec = np.linspace(min(tspan),max(tspan),N)
xtraj = np.zeros((len(x0),N))
row,col = np.shape(xtraj)
impact_ham = []
for i in range(N):
perc = (((i+1)/N)*100)
if i == 0:
print('%6.2f%%' % perc, end='')
else:
print('\r%6.2f%%' % perc, end='')
xtraj[:,i]=integrate_f(f,x,dt,tvec[i])
# check for collision
for p in range(len(list_of_phi_fn)):
if impact_cond(xtraj[:,i],list_of_phi_fn[p]):
H1 = ham_fn(xtraj[:,i-1])
xtraj[:,i] = impact_update_f(xtraj[:,i-1],list_of_p_eq[p])
H2 = ham_fn(xtraj[:,i])
impact_ham.append(np.abs(H2-H1))
x = np.copy(xtraj[:,i])
print(' [COMPLETE]')
return xtraj, impact_ham
def dyn(s,t):
return np.array([*s[int(len(s)/2)::], *[lag_fn[x](s,t) for x in range(6)]])
X0 = np.array([0,0,-np.pi/3,0,0.3,np.pi/6,0,0,0,0,0,0])
tf = 10
dt = 0.01
traj_impact,impact_ham = sim_impact(dyn, X0, [0, tf], dt, integrate, impact_update, impact_condition, phi_fns, p_eq)
100.00% [COMPLETE]
import matplotlib.pyplot as plt
plt.style.use('ggplot')
plt.rcParams.update({'font.size': 14})
fig = plt.figure(figsize=(16, 8))
ax = plt.subplot(211)
ax1 = plt.subplot(212)
ax.plot(np.arange(0,tf,dt),traj_impact[0,:],label='$x_{box}$')
ax.plot(np.arange(0,tf,dt),traj_impact[1,:],label='$y_{box}$')
ax.plot(np.arange(0,tf,dt),traj_impact[3,:],label='$x_{jack}$')
ax.plot(np.arange(0,tf,dt),traj_impact[4,:],label='$y_{jack}$')
ax1.plot(np.arange(0,tf,dt),traj_impact[2,:],label='$\\theta_{box}$')
ax1.plot(np.arange(0,tf,dt),traj_impact[5,:],label='$\\theta_{jack}$')
ax.set_xlabel('time (seconds)')
ax.set_ylabel('distance')
ax1.set_xlabel('time (seconds)')
ax1.set_ylabel('angle (radians)')
ax.set_title('distance over time')
ax1.set_title('$\\theta$ over time')
ax.legend()
ax1.legend()
plt.tight_layout()
plt.show()
fig = plt.figure(figsize=(8, 4))
ax3 = plt.subplot(111)
ax3.scatter(np.arange(1,len(impact_ham)+1,1),impact_ham)
ax3.set_xlabel('Impact Number')
ax3.set_ylabel('|H($\\tau^+$)-H($\\tau^-$)|')
ax3.set_title('Change in Hamiltonian at Each Impact')
plt.tight_layout()
plt.show()
def animate_box_jack(s_array,L=1,l=0.1,T=10):
"""
Function to generate web-based animation of double-pendulum system
Parameters:
================================================
theta_array:
trajectory of theta1 and theta2, should be a NumPy array with
shape of (2,N)
L1:
length of the first leg
L2:
length of the second leg
T:
length/seconds of animation duration
Returns: None
"""
################################
# Imports required for animation.
from plotly.offline import init_notebook_mode, iplot
from IPython.display import display, HTML
import plotly.graph_objects as go
#######################
# Browser configuration.
def configure_plotly_browser_state():
import IPython
display(IPython.core.display.HTML('''
<script src="/static/components/requirejs/require.js"></script>
<script>
requirejs.config({
paths: {
base: '/static/base',
plotly: 'https://cdn.plot.ly/plotly-1.5.1.min.js?noext',
},
});
</script>
'''))
configure_plotly_browser_state()
init_notebook_mode(connected=False)
###############################################
# Getting data from pendulum angle trajectories.
# xx1=L1*np.sin(theta_array[0])
# yy1=-L1*np.cos(theta_array[0])
# xx2=xx1+L2*np.sin(theta_array[0]+theta_array[1])
# yy2=yy1-L2*np.cos(theta_array[0]+theta_array[1])
N = len(s_array[0]) # Need this for specifying length of simulation
# corners of the box
Ba = np.zeros((2,N))
Bb = np.zeros((2,N))
Bc = np.zeros((2,N))
Bd = np.zeros((2,N))
# corners of the jack
J0 = np.zeros((2,N))
J1 = np.zeros((2,N))
J2 = np.zeros((2,N))
J3 = np.zeros((2,N))
for i in range(N):
xb = s_array[0,i]
yb = s_array[1,i]
thb = s_array[2,i]
xj = s_array[3,i]
yj = s_array[4,i]
thj = s_array[5,i]
gwa = g(xb,yb,0)
gab = g(0,0,thb)
gwb = gwa*gab # TRANSFORMATION FOR CENTER OF BOX
gwc = g(xj,yj,0)
gcd = g(0,0,thj)
gwd = gwc*gcd # TRANSFORMATION FOR CENTER OF JACK
t_wb = np.array(gwb).astype(np.float64)
t_wd = np.array(gwd).astype(np.float64)
Ba[:,i] = t_wb.dot(np.array([L/2,L/2,0,1]))[0:2]
Bb[:,i] = t_wb.dot(np.array([-L/2,L/2,0,1]))[0:2]
Bc[:,i] = t_wb.dot(np.array([-L/2,-L/2,0,1]))[0:2]
Bd[:,i] = t_wb.dot(np.array([L/2,-L/2,0,1]))[0:2]
J0[:,i] = t_wd.dot(np.array([l/2,0,0,1]))[0:2]
J1[:,i] = t_wd.dot(np.array([0,l/2,0,1]))[0:2]
J2[:,i] = t_wd.dot(np.array([-l/2,0,0,1]))[0:2]
J3[:,i] = t_wd.dot(np.array([0,-l/2,0,1]))[0:2]
####################################
# Using these to specify axis limits.
xm = -1.5 #np.min(xx1)-0.5
xM = 1.5 #np.max(xx1)+0.5
ym = -1.5 #np.min(yy1)-2.5
yM = 1.5 #np.max(yy1)+1.5
###########################
# Defining data dictionary.
# Trajectories are here.
data=[
# note that except for the trajectory (which you don't need this time),
# you don't need to define entries other than "name". The items defined
# in this list will be related to the items defined in the "frames" list
# later in the same order. Therefore, these entries can be considered as
# labels for the components in each animation frame
dict(name='Box'),
dict(name='Jack'),
dict(name='Jack'),
]
################################
# Preparing simulation layout.
# Title and axis ranges are here.
layout=dict(autosize=False, width=1000, height=750,
xaxis=dict(range=[xm, xM], autorange=False, zeroline=False,dtick=1),
yaxis=dict(range=[ym, yM], autorange=False, zeroline=False,scaleanchor = "x",dtick=1),
title='Jack in Box Simulation',
hovermode='closest',
updatemenus= [{'type': 'buttons',
'buttons': [{'label': 'Play','method': 'animate',
'args': [None, {'frame': {'duration': T, 'redraw': False}}]},
{'args': [[None], {'frame': {'duration': T, 'redraw': False}, 'mode': 'immediate',
'transition': {'duration': 0}}],'label': 'Pause','method': 'animate'}
]
}]
)
########################################
# Defining the frames of the simulation.
# This is what draws the lines from
# joint to joint of the pendulum.
frames=[dict(data=[# first three objects correspond to the arms and two masses,
# same order as in the "data" variable defined above (thus
# they will be labeled in the same order)
dict(x=[Ba[0,k],Bb[0,k],Bc[0,k],Bd[0,k],Ba[0,k]],
y=[Ba[1,k],Bb[1,k],Bc[1,k],Bd[1,k],Ba[1,k]],
mode='lines',
line=dict(color='black', width=2),
),
dict(x=[J0[0,k],J2[0,k]],
y=[J0[1,k],J2[1,k]],
mode='lines',
line=dict(color='blue', width=2),
),
dict(x=[J1[0,k],J3[0,k]],
y=[J1[1,k],J3[1,k]],
mode='lines',
line=dict(color='blue', width=2),
),
]) for k in range(N)]
#######################################
# Putting it all together and plotting.
figure1=dict(data=data, layout=layout, frames=frames)
iplot(figure1)
animate_box_jack(traj_impact,L=1,l=0.2,T=20)
ham_array = []
for i in traj_impact.T:
ham_array.append(ham_fn(i))
plt.rcParams.update({'font.size': 14})
fig = plt.figure(figsize=(10, 6))
ax = plt.subplot(111)
ax.plot(np.arange(0,tf,dt),ham_array)
ax.set_xlabel('time (seconds)')
ax.set_ylabel('Hamiltonian Value')
ax.set_title('Hamiltonian over time')
plt.show()